function [err]=Event_Timer(sf, plot_flag)
% Event_Timer:
%   find, using threshold crossing information, events in the data
%   usually will be applied to GAPFREE data.
%   Plots the following:
%		mean rate vs. time.
%     Interevent interval histogram
%     Joint interevent interval histogram (Isi(n) vs Isi(N+1))

% 
% This returns the modified control worksheet, with a new structure, EVENTS

% Modified 2/10/2000 : changed so that changing rates in the middle of the
% run does not cause subsequent traces to be incorrectly analyzed.
% i.e., all time pointers are recomputed for every trace, depending on the 
% current rate...
% Affected also: find_spikes, spike_sum, etc.
% Derived from CC_IV, STA
%
% 5-Mar-2002
%
% Paul B. Manis, Ph.D.
% pmanis@med.unc.edu
%

global VOLTAGE CURRENT DFILE

global CONTROL

QueMessage('Event_Timer - Starting');

hmp = findobj('Tag', 'Event_Times');
if(ishandle(hmp))
   hc = get(hmp, 'Children');
   delete(hc); % remove the children plots.
end;

dat = [];
time = [];
ivresult = []; % initialize it
do_spike = 0; % assume we want it
nspike = 1; % Which spike to go up to in extraction (1 = first, 2 = first and second... etc)

% various time definitions:
t_pre = 40; % Time to compute average before a spike; normally 25 msec
t_spec_w = 5; % latest time before spike to measure spectrum (e.g., before foot)
t_post = 15;  % window to display AFTER the spike to show ahp
t_delay = 15; % time after last step before the spike before we consider the data on the rising phase valid

min_isi = 2;

[DFILE] = synch_file(DFILE, sf);
if(isempty(DFILE) | isempty(VOLTAGE))
   QueMessage('STA - could not get data');
   return;
end;

% check for units and insert them if they are not present.
if(~check_field(CONTROL(sf), 'I_Unit'))
   CONTROL(sf).I_Unit='pA';
end

if(~check_field(CONTROL(sf), 'V_Unit'))
   CONTROL(sf).V_Unit='mV';
end  

if(~check_field(CONTROL(sf), 'T_Unit'))
   CONTROL(sf).T_Unit='ms';
end
if(~check_field(CONTROL(sf), 'sp_deadwin'))
   CONTROL(sf).sp_deadwin = CONTROL(sf).deadwin;
end
% abstract general information
spike_thresh=number_arg(CONTROL(sf).thresh);
QueMessage(sprintf('Event_Timer - Detection Threshold = %7.2f mV', spike_thresh));
protocol=deblank(lower(CONTROL(sf).protocol));

[records,pts]=size(VOLTAGE);

durs = [CONTROL(sf).durho CONTROL(sf).durs1 CONTROL(sf).durs2 CONTROL(sf).durs3]; 

% compute the time base for plotting (time is [rec, npts] array for EACH record)
time=make_time(DFILE);
tmax=max(max(time));
RATES = DFILE.rate * DFILE.nr_channel / 1000; % array of sampling rates, convert to msec

% get spike times
QueMessage('Event_Timer - Finding Events');
[first_event, first_iei, nr_events, event_train]=find_spikes(DFILE, 0, tmax, spike_thresh);
QueMessage(sprintf('Event_Timer - found %4d spikes', length(event_train.latency)));
first_iei=max(first_iei,0);
%
fsl = first_event(first_event>0);
fisi = first_iei(first_iei>0);
v = size(event_train);
ISILAT=[]; % initialize the ISILAT structure
for i=1:v(2)
   ISILAT(i).nspikes = length(event_train(i).latency);
   ISILAT(i).isi= diff(event_train(i).latency);
   ISILAT(i).lat = event_train(i).latency;
   ISILAT(i).source = event_train(i).source; % save the source trace number too
end

%--------------------------------------------------------------------
% make a figure
%--------------------------------------------------------------------
h = findobj('Tag', 'Event_Timer');
if(isempty(h))
   h = figure('Tag', 'Event_Timer', 'name', 'Event Times/Intervals', 'NumberTitle', 'off');
end;
figure(h);
clf;
set(gcf, 'backingstore', 'off', 'renderermode', 'auto');
fsize = 7;
msize = 1.5;
subplot('Position',[0.0,0.93,1,0.07])
axis([0,1,0,1])
axis('off')
text(0.1, 0.0, sprintf('File: %s   Time: %2d:%2d', DFILE.fullfile, ...
	floor(DFILE.ztime(1)/3600), floor(mod(DFILE.ztime(1), 3600)/60)), 'interpreter', 'none');

subplot('Position', [0.10, 0.20, 0.6, 0.75]);
% split data into multiple traces on the plot...
maxt=15000; % 15 second stretches...
np = ceil(max(time)/maxt); % find number of plots needed.
voff = 0; % do it with a voltage offset...
minv = min(VOLTAGE);
maxv = max(VOLTAGE);
dt = time(2)-time(1); % assume constant sampling.
nt = floor(maxt/dt);
voff = 1.05*abs(maxv-minv); % find range so no overlap is present
	p = 1/length(v(2));
	st = spike_thresh*ones(length(event_train.latency),1);
for i = 1:np
   t0 = (i-1)*maxt;
   t1 = t0+maxt;
   p_t0 = (i-1)*nt+1;
   p_t1 = p_t0 + nt-1;
   plot((time(p_t0:p_t1)-t0)/1000, VOLTAGE(1,p_t0:p_t1)'-(i-1)*voff, '-k');
	hold on;
   sti = find(event_train.latency >= t0 & event_train.latency < t1);
   % mark where spikes detected
	plot((event_train.latency(sti)-t0)/1000, st(sti)-(i-1)*voff, 'linestyle', 'none', 'marker', 'o', ...
      'color', 'r', 'markerfacecolor', 'none', 'markersize', msize);
end;
xlabel('time (sec)');
   ylabel('V (mV)');
   
subplot('Position', [0.1, 0.1, 0.6, 0.07]);
for i = 1:v(2)
   x = [ISILAT(i).lat];
   y = [ISILAT(i).isi];
   plot(x(1:end-1)/1000, y(1:end)/1000, 'linestyle', 'none', 'marker', 'o', ...
         'color', 'k', 'markerfacecolor', 'black', 'markersize', msize);
end;
xlabel('time (sec)');
yLabel('ISI (sec)');

subplot('Position', [0.75, 0.65, 0.23, 0.23]);
hist([ISILAT.isi]/1000, 20);
xlabel('ISI (sec)');
ylabel('N');
title('ISI distribution')
subplot('Position', [0.75, 0.1, 0.23, 0.23]);
plot([ISILAT.isi(1:end-1)]/1000, [ISILAT.isi(2:end)]/1000, 'linestyle', 'none', 'marker', 'o', ...
   'color', 'k', 'markerfacecolor', 'black', 'markersize', msize);
u1 = get(gca, 'YLim');
u2 = get(gca, 'XLim');
u1(1) = 0;
u2(1) = 0;
x = max(u1(2), u2(2));
u1(2) = x; u2(2) = x;
set(gca, 'YLim', u1);
set(gca, 'XLim', u2);
hold on;
plot(u1, u2, 'color', 'k');
title('Joint ISI distribution');
xlabel('ISI1 (sec)');
ylabel('ISI2 (sec)');
return;
